Leçon 8: Différences finies avec numpy.

Exemple de solution de transport diffusion-convection scalaire, la vitesse du fluide est connue, on ne solutionne pas Navier-Stokes.

On utilise la solution matricielle directe des équations de différences finies.

L'équation est donc $$ v_x\frac {\partial C}{\partial x } +v_y \frac {\partial C}{\partial x }= D \left (\frac {\partial ^{2} C}{\partial x^{2} } +\frac {\partial ^{2} C}{\partial y^{2} } \right)$$

On utilise des différences décentrées arrières pour la convection, donc les vitesses données doivent être positives. On ne discutera pas ici des problèmes associés au décentrage et à la stabilité des équations de diffusion-convection. Il serait toutefois assez facile d'appliquer un algorithme de décentrage de type 'upwind' qui assurerait la stabilité (sinon la précision) des équations.

L'exemple donné est adapté de Heys, Chemical and Biotechnological Engineerinc Calculations using Python 3.

#
# Une fonction pour construire la matrice de D.F.
#
def differencesFinies(NPI,NPJ,aid,ajd,aic,ajc,ap,A,b):
    for i in range( 0 ,NPI):
        for j in range( 0 ,NPJ) :
            p = j*    (NPI)    + i
            pE= j*    (NPI)    + i+1           #  Construction de la matrice de différences finies A x = b
            pO= j*    (NPI)    + i-1           #  de taille (NPI x NPJ)  par (NPI x NPJ)
            pN= (j+1)*(NPI)    + i 
            pS= (j-1)*(NPI)    + i 
            if   i == 0:
                A[p,p] =1.0
                b[p]= np.sin(np.pi*y[j])
            elif i == NII:
                A[p,p]=1.0                      # conditions aux limites
                b[p]=0.0
            elif j == 0: 
                A[p,p]=1.0
                b[p]=0.0
            elif j == NIJ: 
                A[p,p]=1.0
                b[p]=0.0
            else:
                A[p,p]  = ap+aic+ajc             # matrice, points intérieurs
                A[p,pE] = aid
                A[p,pO] = aid-aic
                A[p,pN] = ajd
                A[p,pS] = ajd-ajc
                b[p]    = 0
    return A,b 
#
#
# Début du code
#
import numpy as np
from numpy.linalg import solve
import pylab
%matplotlib inline
NII     = 15        # nombre d'intervalles direction x
NIJ     = 20        # nombre d'intervalles direction y
NPI     =NII+1      # nombre de points direction x
NPJ     =NIJ+1      # nombre de points direction y
ouest = 0.0 
est   = 1.0
sud   = 0.0
nord  = 1.0
x     = np.linspace(ouest , est  ,NPI )
y     = np.linspace( sud  , nord ,NPJ )
v     = np.array([5,5])
dx    = x[1]-x[0]
dy    = y[1]-y[0]
aic   = -v[0]/dx
ajc   = -v[1]/dy
aid   = 1.0/dx**2
ajd   = 1.0/dy**2
ap   = -2.0*(aid+ajd)
X,Y   = np.meshgrid(x,y)
total = (NPI)*(NPJ)
#
#     initialisation du système
#
A     = np.zeros((total,total))
b     = np.zeros(total)
#
#     
A,b   = differencesFinies(NPI,NPJ,aid,ajd,aic,ajc,ap,A,b)
#
#
c=solve(A, b )                   #   solution directe
C=c.reshape(NPJ,NPI)             #   reconstruire la matrice pour tracer
pylab.contourf(X,Y,C)
pylab.colorbar ( )
pylab.xlabel ('x')
pylab.ylabel ('y')
Ti=' Diffusion-convection de C'
pylab.title(Ti)
pylab.show( )